Lyapunov indices with two nearby trajectories in a curved spacetime 
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We compare three methods for computing invariant Lyapunov exponents (LEs) in general rela- 
tivity. They involve the geodesic deviation vector technique (Ml), the two-nearby-orbits method 
with projection operations and with coordinate time as an independent variable (M2), and the 
two-nearby-orbits method without projection operations and with proper time as an independent 
variable (M3). An analysis indicates that Ml and M3 do not need any projection operation. In 
general, the values of LEs from the three methods are almost the same. As an advantage, M3 is 
simpler to use than M2. In addition, we propose to construct the invariant fast Lyapunov indictor 
(FLI) with two-nearby-trajectories and give its algorithm in order to quickly distinguish chaos from 
order. Taking a static axisymmetric spacetime as a physical model, we apply the invariant FLIs 
to explore the global dynamics of phase space of the system where regions of chaos and order are 
clearly identified. 
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I. INTRODUCTION 

Chaos is often visible in nonlinear systems. General 
relativity as a nonlinear theory is potentially chaotic. Al- 
though many features of chaos in Newtonian dynamics 
have been known clearly for over forty years, their ap- 
plication in relativistic astrophysics began to be widely 
appreciated only within the last decade or so [1]. One 
main interest lies in studying the difference of dynamics 
between Newtonian and relativistic trajectories. Varadi 
et at [2] noted that the general relativity effects are 
small for the outer planets but not negligible. Addition- 
ally, Wanex [3] revealed the chaotic amplification effect 
in the relativistic restricted three-body problem (namely, 
the ideal spacecraft-earth- moon orbital system ) . In par- 
ticular, several authors found chaos in two relativistic 
systems including the two fixed black holes [4] and the 
Schwarzschild black hole plus a dipolar shell [5,6], which 
does not appear in their corresponding Newtonian coun- 
terparts at all. Another important problem that has 
been developed in recent years is to question whether 
spinning compact binaries [7-15] as promising sources of 
gravitational waves exhibit chaotic behavior because the 
gravitational- wave detection can not succeed when chaos 
is present. In practice, all the examples are attributed 
to the geodesic or non-geodesic motion of particles in 
a given gravitational field. On the other hand, the time 
evolution of the gravitational field itself, such as the mix- 
master cosmology [4,16,17], is also of great interest. 

In general, the methods for quantifying the ordered 
or chaotic nature of orbits in general relativity follow 
those widely applied in classical physics. Let us recall 
a fraction of these classical methods briefly. Poincare's 
surface of section is one of the most common qualitative 
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tools in the analysis of conservative Newtonian dynam- 
ical systems of not more than 2 degrees of freedom or 
with 3 degrees of freedom and axial symmetry. However, 
this technique is difficult to describe a higher dimensional 
phase space. Certainly, the principal Lyapunov exponent 
(LE), as a measure of the average exponential deviation 
of two nearby orbits, is frequently used. There are two 
different algorithms for the calculation of LE: the vari- 
ational method and the two-particle one [18] (see Sec. 
II) . The algorithm of LE is applicable to a phase space 
with any dimension, but a high dimension would cause 
extremely expensive computation to get a reliable value 
of LE. There are also other qualitative methods for multi- 
dimensional systems, such as the power spectra, smaller 
alignment index (SALI) [19,20] and fast Lyapunov indica- 
tors (FLIs) [21,22] etc. The power spectra display a finite 
number of discrete frequencies for regular orbits, whereas 
they are continuous for chaotic ones. One of its problems 
is that it is ambiguous to differentiate among complicated 
periodic orbits, quasi-periodic orbits and weakly chaotic 
orbits. As far as the SALI method is concerned, it is 
based on the difference or the sum of two normalized de- 
viation vectors at the same points of an orbit. If the 
dimension of phase space is larger than 2, the indicator 
changes around non-zero values for regular orbits, while 
it tends exponentially to zero for chaotic orbits. The 
FLI uses the logarithm of the length of a deviation vec- 
tor, which increases following completely different time 
rates for different orbits thus allowing one to distinguish 
between the ordered and chaotic cases (see Sec. Ill for 
some details). Both the SALI and the FLI are not only 
very fast tools to find chaos, but also can sketch out the 
global structure of phase space. Obviously it takes more 
CPU time for the former than for the latter since the for- 
mer requires two deviation vectors, while the latter needs 
only one. 

Considering the above various cases, we are mainly 
in favor of applying two Lyapunov indices, the LE and 
the FLI, to study the geodesic or non-geodesic motion 
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of particles in relativistic gravitational systems with at 
least two degrees of freedom. As is well known, the clas- 
sical definition of LE depends on the choice of time and 
space coordinates in general relativity, because a coor- 
dinate gauge in relativity can be arbitrarily adopted so 
that time and space coordinates are not necessarily phys- 
ical. In other words, there would be different values of 
the classical LEs in different coordinate systems. Even 
such LEs may provide wrong information on the dynam- 
ical features of a system. For instance, the maximum LE 
of a chaotic system turns out to be zero after a logarith- 
mic time transformation, that is to say, chaos becomes 
hidden. Consequently, it is strongly desirable to develop 
a coordinate independent definition of LEs for the study 
of relativistic dynamics. 

In general relativity, one can still follow the two ap- 
proaches to calculate LEs as in the Newtonian case. As a 
counterpart of the variational technique, a geodesic devi- 
ation vector can be obtained from the geodesic deviation 
equation for a given geodesic flow in general relativity. 
Using it, one can easily get a method (Ml) for the in- 
variant LEs in the configuration space as well as in the 
phase space. In addition, Imponente and Montani [17] 
presented an invariant treatment by projecting a geodesic 
deviation vector for the Jacobi metric on an orthogonal 
tctradic basis. In this way they can succeed in gaining 
an insight into the dynamics of the mixmaster cosmology. 
On the contrary, it is seldom to see in the references that 
LEs are directly computed by use of the geodesic devi- 
ation equation. Perhaps it is rather troublesome to de- 
rive the complicated curvature tensor. In this sense, it is 
very natural to extend the two-particle method from the 
classical to the relativistic case. It is convenient to em- 
ploy the method (M2) introduced by Wu and Huang [23]. 
They compute gauge invariant Lyapunov exponents by 
calculating the space separation between an "observer" 
and a "neighbor" particles that move along two neigh- 
boring orbits in the phase space. Here it is necessary 
to make use of a "1+3" split of the observer's spacetime 
and its space projected operator. Meanwhile, Wu and 
Huang [23] use coordinate time as the time variable of 
the equations of motion. Besides M2, there is another 
invariant two-particle method (M3), in which the inte- 
gration time variables are the individual proper times 
of the two particles instead of the common coordinate 
time and the projection operation is not used. A method 
like M3 can be seen in many references (for example, see 
[8,12]) where the proper time and the Euclidian distance 
in the phase space are adopted, but we prefer to utilize 
the Riemannian distance in the configuration space than 
the Euclidian one in the phase space for conceptual clar- 
ity. 

Tancredi et al. [24] have compared the two approaches 
for the calculation of LEs in Newtonian dynamics. It 
is still necessary to do a comparison of the three meth- 
ods because general relativity is much different from the 
Newtonian dynamics. A fundamental task is to select 
an optimal and valid algorithm to compute LEs. An- 



other aim of our work is to provide a simple, rapid and 
efficient indicator of chaos which is independent of the di- 
mension of phase space and the choice of time and space 
coordinates. It is assumed that this indicator not only 
can tell us some information about the global motion in 
a complicated multi-dimensional relativistic system, but 
also study the transition from regular motion to chaos as 
certain physical parameters alter. 

In principle, it is possible to fulfill the two purposes in 
terms of LE. Udry and Pfenniger [25] made a detailed 
quantitative estimate of chaos in a series of 24 triax- 
ial models of elliptical galaxies by means of LEs. They 
claimed that an interval of about 2 Hubble times corre- 
sponding to 66 periods or so is sufficient to get stable 
values of LEs for each of one hundred randomly selected 
orbits in every model. However, as Contopoulos and Bar- 
banis [26] found, lots of these orbits must be calculated 
for 12500 periods, and sometimes up to 250000 periods 
for the accomplishment of reliable values of LEs. On the 
basis of the values of the LEs, they displayed the struc- 
ture of phase space of the systems. In addition, Carani- 
colas and Papadopoulos [27] studied the transition from 
regular motion to chaos in a 2-dimcnsional logarithmic 
potential of elliptical galaxies by observing the Poincare 
surface of section as some dynamical parameters change. 
Of course, it is necessary to calculate thousands of orbits 
in the above cases, and each of them must be computed 
for long enough time in general. Consequently, the time 
required to compute such numbers is rather expensive, 
even would reach the limit of the present numerical ex- 
periments. Under these circumstances, one had better 
take advantage of the FLI as a quicker and more sensitive 
indicator to separate chaotic orbits from regular ones. 
As is mentioned above, the indicator is originally based 
on the tangential vectors from the variational equations 
of Newtonian dynamics [21]. Pondering the difficulty 
and complicacy to derive variational equations in general 
relativity, we shall be interested in developing the FLI 
method with the two particles approach — the FLI com- 
puted with two-nearby-orbits. It should be emphasized 
that this approach would meet difficulty without renor- 
malization in the numerical calculation, different from 
the original approach of FLIs [21] with the tangential 
vector. An algorithm of FLI in the present paper will be 
provided. 

This paper is organized as follows. The classical def- 
inition of LE is reviewed briefly, and the three methods 
Ml, M2 and M3 in relativistic systems are introduced in 
Sect. II. Meanwhile, we offer a theoretical analysis of Ml 
and M3 to explain why they do not need projection oper- 
ations. In Sect. Ill the FLI with two-nearby-orbits and 
its algorithm are presented. Setting a static axisymmet- 
ric spacetime composed of a Schwarzchild black hole and 
an octopolar shell as a test model, we investigate the 
global dynamics of this system. Finally, the summary 
follows in Sec. IV. 

Throughout the work we use units c = G = 1, and 
take the signature of the metric as (—,+,+,+). Greek 
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subscripts run from to 3 and Latin indexes from 1 to 
3. The symbol o denotes the Euclidian inner product; 
while \a\ is the length of the vector a in the Euclidian 
space. The symbol \a\ represents the absolute value of 
a scalar a. In addition, • stands for the Riemannian 
inner product, and ||£|| is the Riemannian norm of the 
vector £ corresponding to the tensor £ Q or £ Q . / is a set 
of functions. We specify t as a time coordinate variable, 
and t, a proper time variable. D/Dt denotes a covariant 
derivative operator vs t. 



II. COMPARISONS OF THREE METHODS FOR 
COMPUTING LE IN GENERAL RELATIVITY 

A. LE in classical physics 

In classical physics the LE is to characterize the mean 
exponential rate of divergence of trajectories surrounding 
a given trajectory in the phase space. For a compact 
autonomous n-dimensional system 



f(X) 



(1) 



with the solution X(t) = (x,x) and its corresponding 
variational equation 



(2) 



with a solution Y(t) 
is given by 



(£, £) at time t, the maximum LE 



7l =hm 



(3) 



This is called as the variational method. The system is 
chaotic if 71 > 0, otherwise it is regular. This way is 
rigorous to get the tangent vector Y(t), but it is rather 
cumbersome to derive Eq. (2) in general. Because of that 
people usually use the deviation vector AX(t) between a 
reference trajectory X{t) and a shadow one X(t) as the 
approximate tangent vector. This is a less rigorous but 
still useful technique, so-called the two-particle method 
[18,24] or two-nearby-trajectories method, in which the 
LE in the expression of Eq. (3) is to be replaced by 



.. 1. \AX(t)\ 
72 = hm - In — — — - . 

' t^oo t I AX (0)| 



(4) 



Here an initial separation |AX(0)| relative to |X(0)| not 
larger than 10 -8 is viewed as the best choice to guarantee 
AX as a good approximation to Y and to avoid the 
overestimation of LE [24] . 

As was mentioned in Ref. [23] , it is preferable to com- 
pute the LE in the configuration space instead of in the 
phase space as LEs in the two spaces are both effective 
in detecting the long-term dynamical behavior of orbits. 



In this case, Eqs. (3) and (4) are, respectively, modified 
as follows 



Ai 
A 2 



|€(o)|' 

t^oo t |Aaj(0)| 



(5) 
(6) 



For the two particles method, it is necessary to scale 
the distance |Aa:(f)| down from time to time. In this 
way, the shadow trajectory returns to the neighborhood 
of the reference one along the deviation vector Ax(t). 
The magnitude of Ax(t) shrinks to the initial distance 
|Ax(0)| after each renormalization, and the velocity de- 
viation vector Ax(t) should be multiplied by the same 
factor |Aa;(0)|/|Acc(t)|. Meanwhile, it is also vital to 
avoid saturation of orbits in a bounded chaotic region. 

The time t and coordinate x in Eq. (5) or Eq. (6) are 
not necessary physical and meaningful in general relativ- 
ity. Thus, it is desirable to define covariant LE. 



B. LEs in general relativity 



For a given 4-dimensional spacetime with the metric 
ds 2 = g^ v dx^dx u , first we consider a rigorous definition 
of LEs by using a geodesic deviation vector. 



1. Geodesic deviation vector technique 

In the spacetime, a free particle moves along the 
geodesic equation DU^/Dt = 0, that is to say, 



^ +<* -p 



and its geodesic deviation equation is 



(7) 



(8) 



where T^, U a ( 



stand for the Christoffel 



x°) and 

symbol, 4- velocity, and the Riemannian curvature tensor, 
respectively. 

From Eq. (5), it is easy to get an invariant definition of 
LE if the Riemannian norm and proper time r substitute 
the Euclidian norm and coordinate time t, respectively. 
In the geodesic deviation vector method (Ml) the LE is 
given in the form 



Ai = IhnilnMM 
r^cor ||£(0)| 



Uir) 



(9) 
(10) 



As an illustration, a similar method given in the phase 
space can be seen in Ref. [28]. It is easily shown that 
the Riemannian inner product £ • £ is positive definite 
because £(t) is always space-like when £(0) is space-like 
for a given geodesic flow. 
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On the other hand, there is not a variational equation 
like the form (8) for a non-geodesic flow. However, it is 
still possible to derive the variational equation similar to 
Eq. (2) for the non-geodesic flow. For example, Nieto 
et al. [29] gave a relativistic top deviation equation as 
a generalization of the geodesic deviation equation for a 
pair of nearby point particles. Thus Ml remains valid if 
the spacetime background is ds 2 = g^dx^dx" . 



2. Two-nearby-orbits method with projection operations 

The derivation of the geodesic deviation equation is 
usually a rather hard task. In particular, it is much ardu- 
ous to obtain the variational equation of a non-geodesic 
flow. Thus, it is a good idea to refine the classical defi- 
nition by Eq. (6) using an invariant version of it, intro- 
duced by Wu and Huang [23] (method M2). 

Let two particles, an "observer" and his "neighbor", 
move on two nearby trajectories in a curved spacetime. 
At a coordinate time t the observer is at the point O 
with coordinate x a and 4- velocity U a , and his neighbor 
reaches the point O with coordinate x a . The deviation 
vector 



Ax(t) = Ax a (t) = x a (t) - x a (t) 



(11) 



from O to O is projected to the observer's local space and 
the resulting projected vector is Ax" = /iSAar, where 
h a/3 = g a " + U a UP is the space projection operator of the 
observer. The space distance of the neighbor measured 
by the observer at time t is 



A7 = \/ g a0 AxlAx>l = yfh aP Ax a AxP. (12) 



Hence we define an invariant LE (M2) as 



. ,. 1, AL(r) 
A 2 = hm -In- — 

t^oo T AL(0) 



(13) 



where the proper time t corresponds to the coordinate 
time t according to the metric. Obviously, this LE is 
independent of coordinate transformations [23]. 

There are detailed implementations of M2 in [23] . Here 
we emphasize that the coordinate time is adopted as the 
common independent variable for both particles in their 
equations of motion and one has to construct an equation 
for dr/dt, which is to be integrated together with the 
motion equations to get the proper time of the observer. 
The reason for doing this is that the two particles have 
their own and different proper times but the numerical 
integration demands a common time variable. 



conditions with the proper time as an integration vari- 
able, we attain the deviation vector 



Ax(t) = Ax a (t ) = x a (r) - x a {r) 



(14) 



at the proper time t of the observer. Here a difference 
between Eq. (11) and Eq. (14) should not be neglected. 
The difference is that Ax is a function of t in Eq. (11), 
while being a function of r in Eq. (14). In this case, 
we can give another two-nearby-orbits method (M3) as 
follows 



r^oo T \\Ax 



(15) 



||Asb(t)|| = ^\Ax.Ax\ = yt^l. (16) 

Next let us study the three methods from the theoret- 
ical and numerical points of view. 



C. The reason why Ml and M3 do not need 
projection operations 

Obviously, projection operations do not appear in Ml. 
Similarly, Hartl noted this fact in his PhD thesis [8]. 
Now, we give a discussion on this issue. It has been shown 
without any difficulty that £(t) is always perpendicular 
to the 4- velocity U(t) at any (proper) time provided that 
£(0) and D£(0)/Dt are both perpendicular to U(0), re- 
spectively [28]. What would happen if £(0) • 17(0) ^ 
or D£(0) I Dt • 17(0) ^ 0? Let us explore the question in 
the following. 

Noting that £ • U and D£/Dt • U are scalars, we can 
operate a covariant derivative in place of an ordinary one. 
The demonstration is 



dv Dt ' 



Dt ' Dt ? Dt 



D 2 £ 



Rv a xpU a U l3 U»e- 



(17) 



(18) 



The result of Eq. (18) is obviously identical to zero be- 
cause it is both symmetric and antisymmetric with re- 
spect to the indicies v and a. Consequently, it can be 
inferred that the right side of Eq. (17) is a constant C, 
which turns out to be the value of D£(Q)/Dt* 17(0). As 
was stated above, £»J7 = 0ifC = 0. Otherwise we have 



£ • U = Ct + C, 



(19) 



3. Two-nearby-orbits method without projection operations 

On the other hand, if we integrate directly the equa- 
tions of motion, Eq. (7), for two slightly distinct initial 



where C is another constant corresponding to the initial 
value of £ • U. Now let us investigate the physical sig- 
nificance of |£ • U\. Obviously it is no other than the 
length of the projected vector — (£• U)U, the time com- 
ponent of £ as measured by the particle. Equation (19) 
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shows that the evolution of £ along the time direction 
of the particle is linear, therefore, the LE in this direc- 
tion vanishes. In fact, the result is very natural due to 
the integral U a U a = —1. Therefore, £ and its projected 
vector along the space direction of the particle change at 
the same rate with its proper time. This implies that it 
is not necessary to project £ into the space direction in 
the calculation of the principal LE. In other words, it is 
reasonable to use Ml as the standard definition of the 
maximum LE in general relativity. 

As to M3, the deviation vector Ax(t) in Eq. (14) 
is viewed as an approximation to the geodesic deviation 
vector £ by Eq. (8). Therefore, M3 is nearly the same as 
Ml as long as the length of £ is kept small enough. 

For M2, we use the projected vector As" rather than 
Ax(t) used in Eq. (11) as an approximation of Ax(t) 
in Eq. (11) is the deviation vector between the observer 
and his neighbor at the same coordinate time t, therefore 
it is necessary for Ax(t) to be projected along the space 
direction of the observer. On the other hand, Ax(t) and 
U do not necessarily obey the relation (19) when t and 
r have large difference. Then in the observer's in-track 
direction he may find LE in the form 

1 AT (t^) 
\n-track = lim - In , (20) 

T^OO T Al (0) 

where AT(t) = \Ax(t) • U\ is the projection of Ax(t) 
along the observer's time direction. This shows that one 
may find chaos in the observer's time direction if \n-track 
is used as an index for chaos. This displays it necessary 
to employ the projection norm for M2. 

We shall further check the validity of the three methods 
through practical calculations. 

D. Numerical comparisons 

Let us take a Weyl spacetime as a physical model (see 
the appendix), and choose parameters as E = 0.9679, 
L = 3.8 and O = 7.012 x 10~ 7 . The initial conditions of 
two variables r and r are arbitrarily given within their 
respective admissible intervals except 9 = tt/2, while 6 is 
derived from Eq. (A4). First order differential systems 
from the geodesic Eqs. (A5) and (A6) are computed with 
a 7-8 order Runge-Kutta-Fehlberg algorithm of variable 
time-step (RKF7(8)). A proper time output interval is 
0.1, and the Poincare surface of section is at the plane 
9 = ^ (9 < 0) in Fig. 1. The intersections of the orbits 
by this Poincare surface of section describe the global 
dynamical feature of the system. There are two regions 
in the (r, f) plane. One region with randomly distributed 
points is regarded as a chaotic one, while the other, a 
regular one consisting of many tori and islands. 

Now we quantitatively describe the dynamics of par- 
ticles moving in this spacetime by employing the LEs 
above. For Ml, it is necessary to integrate the dynam- 
ical Eqs. (A2)-(A6) numerically with their geodesic de- 
viation equations together so that we obtain their solu- 



tion (x a , U a ) and their variational solution (£,£) in the 
forms | = (St,8r,89,5(j)) and £ = (Si,Sr,S9,5<p). Here 
the proper time r is chosen as an integration variable. 
In addition, the renormalization technique is used if nec- 
essary. As far as M2 is concerned, we numerically trace 
the trajectories of the observer and his neighbor with 
the Eqs. (A7,A8) and the equations involving dr/dt and 
d<fi/dt with slightly different initial conditions. The ini- 
tial conditions of the observer are the same as the above. 
As to his neighbor, an initial separation in the order of 
10~ 8 is adopted only on the r direction, regarded as the 
best choice [24], and the other coordinates remain the 
same as the observer's except 9. This process also fits 
M3, but only Eqs. (A2)-(A6) are integrated numerically. 



Choosing many ordered and chaotic orbits as numeri- 
cal tests, we found that the three methods Ml, M2 and 
M3 almost give the same final values of LEs for each 
orbit. There are only a few differences on the tran- 
sition phases of the orbits. For example, we choose a 
chaotic orbit with the initial conditions: t = 0, 9 = tt/2, 
r = 9, 4> — 0, r — 0.025 and the initial value of 9 derived 
from Eq. (A4). Let the initial deviation vector satisfy 
5t = 66 = 5<f> = Sr = 56 = and 8r = 1/y/gU. As 
is shown in Fig. 2, M3 is always close to Ml (that is, 
M3 is just a good approximation of Ml), but M2 is not 
consistent with Ml until the stabilizing value of LEs is 
obtained. Consequently, the three methods, Ml, M2 and 
M3, have almost the same final value of LEs. 



For conceptual clarity it is of the physical significance 
to define LE as a coordinate gauge invariant in general 
relativity. For this purpose, there are the above three 
methods (Ml, M2 and M3) to compute LEs. Ml is more 
rigorous than M2 or M3. On the contrary, M2 is more 
convenient to use than Ml, and M3 is the simplest. In 
addition, M2 and M3 are easier to treat a non-geodesic 
flow as well as a geodesic flow than Ml. Numerical exper- 
iments display that Ml, M2 and M3 not only can obtain 
the same final values but also have small differences in 
the transition phase for the calculation of LEs. Consid- 
ering the application of convenience, we conclude that 
M3 should be worth recommending to calculate LEs in 
a relativistic gravitational system. It is useful to accept 
some advice on a suitable choice of the initial separation 
and a renormalization time step introduced by Tancredi 
et al [24]. 



Although the LEs have been widely used to distinguish 
between regular and chaotic orbits, their stabilizing val- 
ues obtained often need rather expensive computational 
cost in exploring the global structure of phase space. In 
view of this, we shall modify M3 and adopt its corre- 
sponding FLI in the following section. 
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III. DESCRIPTION OF THE STRUCTURE OF 
PHASE SPACE BY USING FLI WITH 
TWO-NEARBY-TRAJECTORIES 

A. Fast Lyapunov indicators 

The time interval that is necessary to reach a given 
value of either the length of a tangential vector or the 
angle between two tangential vectors can be taken as an 
indicator of stochasticity of a Newtonian dynamical sys- 
tem. Following this idea, Froeschle et al. [21] defined 
three different FLIs in 1997. However two of these meth- 
ods need to solve the variational Eq. (2) for n times when 
a set of n independent initial tangential vectors with the 
same initial conditions are chosen. In this circumstance, 
Froeschle and Lega [22] improved the FLI as follows 



m = ^\Y(t)\, 



(21) 



where Y is a tangential vector from Eq. (2) and | V(0)| = 
1. Given a threshold, the indicator ip reaches a value fast 
for a chaotic orbit, but it would take a rather long time 
for an ordered one. Conversely, in the same time interval 
the indicator tends to different values for ordered and 
chaotic orbits. Namely it grows exponentially with time 
for the chaotic orbit, but only algebraically with time in 
the regular case. This allows one to distinguish between 
the two cases. There is a close relation between the FLI 
(tp) and the LE [Eq. (3)]: the FLI divided by time t tends 
to the LE when the time is sufficiently large. Besides, 
overflow of the lengths of tangential vectors in the case 
of a chaotic orbit can be avoided because the integration 
time is not long enough. This is the reason why the 
indicator is classified as a "fast" method. 

Following the idea above and the coordinate gauge in- 
variance, in a curved spacetime we could easily define an 
invariant FLI corresponding to the LE in Eqs. (9) and 
(15). As far as the geodesic deviation vector method Ml 
is concerned, its corresponding FLI is 



= log 10 U(r) 



(22) 



In the light of Eq. (15), we define the FLI with two- 
nearby-trajectories as follows: 



FLI(t) = log 10 



|As(t)| 
|A*(0)| 



(23) 



Next we shall mainly describe the numerical implemen- 
tation of the FLI and test its validity in the description 
of dynamics. 



B. The algorithm in detail 

The original FLI proposed by Froeschle and Lega [22] 
requires to compute the expansion rate of a tangential 
vector in the variational equations and it does not need 
any renormalization. We now reform it in general rela- 
tivity into a coordinate invariant and compute it with the 



two particles approach. The FLI(t) in Eq. (23) can not 
be computed without renormalization because the dis- 
tance between the two particles, ||Ax||, would expand so 
fast in the case of chaos as to it could reach the chaotic 
boundary to cause saturation. 

In our numerical model (see the appendix) we choose 
||Ax(0)|| = 1(T 9 . We found saturation when ||Ajc|| = 1, 
therefore we choose ||Ax|| = 0.1 as the critical value to 
implement the renormalization. In this way, the number 
of renormalization for computing FLI is less than that 
for LE. This brings an advantage to guarantee the speed 
of computation. Let k (k = 0, 1, 2, • • • ) be the sequential 
number of renormalization, then we calculate the FLI 
with the following expression 



FLI k (r) - -fc-(l + logic l|A:z(0)||) 
|Ax(r)|| 



+ logi 



|Aas(0)|| 



(24) 



where ||Ajc(0)|| < ||Ajc(t)|| < 0.1. This technique de- 
pends on the choice of the initial deviation ||Ax(0)||, 
and fails when ||Aa;(0)|| is too small or too large. After 
many numerical experiments, we found that ||Aa;(0)|| s=s 
10~ 7 - 10~ 9 works well. 



C. Numerical tests of FLI 

We take the initial separation Ar = 10~ 9 and choose 
the regular orbit M and the chaotic one N in Fig. 1 
to test the sensitivity of the FLI as well as ^. As is 
expected, Fig. 3 (colored light gray) displays that there 
is a drastic difference of the FLIs between the regular 
orbit and the chaotic one. For the ordered orbit, within 
a time span of 10 5 (equivalent to about 100 periods), the 
FLI is smaller than one. However, we clearly see a sharp 
increase up to 24 of the FLI for the chaotic orbit. In 
other words, the length of the geodesic deviation vector 
£ reaches the value 10 24 . As mentioned above, the vector 
increases in an algebraic law for quasi-periodic orbits, but 
it does with an exponential law for chaotic orbits. Indeed 
the chaos of the orbit in the right of Fig. 3 becomes 
explicit after a time span 4.5 x 10 4 . The FLI is a cheaper 
way to distinguish between chaotic and regular orbits 
and explore the global qualitative structure of the phase 
space of a system. 

Another point to emphasize is that the values of FLI 
are very analogous to that of "J/ (see black dots in Fig. 3) . 
The values of FLIs by use of the two methods for the reg- 
ular case have no difference in a long time since renormal- 
ization is not used. For the chaotic case we use only two 
times of renormalization. In general, it is enough to take 
small numbers of renormalization to calculate the FLIs 
in the orbits with weak chaos. This shows sufficiently the 
validity of our FLI with two-nearby-trajectories. 

To examine the validity of the FLI further, in the left 
of Fig. 4 we have plotted the running maximum of FLIs 
as a function of the initial action r for a set of 1657 or- 
bits that are regularly spaced in the interval [5.77, 22.34] 
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on the axis r — in Fig. 1. Here each orbit is inte- 
grated until the time reaches 10 5 . As is described in 
Fig. 1, the regular region is located in one interval be- 
tween i?(8.45,0) and 5(21,0) and another interval be- 
tween #(22.15, 0) and the point (22.34, 0) on the straight 
line r = 0. Obviously we find that all the FLIs in the two 
intervals are not larger than 6, but all the FLIs in the two 
intervals [5.77, 8.45] and [21, 22.14] on the straight line 
f = are larger. This method displays the set of chaos, 
although it does not give the final values of LEs. In a sim- 
ilar way, the orbits along the straight line r = 9 in Fig. 1 
are also separated into regular and chaotic intervals (sec 
the right of Fig. 4). These facts show sufficiently that 
the FLI is rather successful to distinguish between the 
regularity and the chaoticity. However, we must point 
out a problem that the FLI does not provide more de- 
tails about ordered orbits containing quasi-periodic and 
resonant periodic orbits. Here we give an illustration in 
detail. 

On the basis of Froeschle and Lega [22] , the orbit with 
the starting point M(15, 0) in Fig. 1 corresponds to a res- 
onance since the maximum value of FLI arrives at about 

0. 68 as the smallest value in the left of Fig. 4. However 
this is not a resonant periodic orbit but just a quasi- 
periodic one. Actually the point 7(14.00349, -0.0312004) 
of the Poincare map in the middle of tori is a fixed 
point. We found that all tori between the torus M 
and the point I have almost the same values of FLI. 
We cannot say that FLI is appropriate to identify the 
resonance categories. As another example, let us fo- 
cus a trajectory made up of three little loops or islands 
surrounding three invariant points J\ (9.0975, —0.0281), 
J 2 (17.13931, 0.04913) and J 3 (13.94246, -0.08936) in Fig. 

1, respectively. An important point to mention lies in the 
manner where the three little loops appear on the plot as 
the same trajectory. Rather than tracing one loop at a 
time, successive points occur at each of the three loops in 
turn. In this sense, the three fixed points are inhabited in 
the same periodic orbit. The maximum value of FLI we 
computed for this periodic (resonant) orbit after a time 
span of 10 5 turns out to be 1.75 or so. Obviously the 
value of FLI for the resonant orbit is larger than that for 
the quasi-periodic orbit M. In a word, it is difficult to 
apply values of FLI to distinguish various regular orbits. 



D. Classification of orbits by FLI 

Many experiments above have shown that our FLI is a 
very sensitive tool for detecting the regular and chaotic 
orbits. Next we shall follow this FLI to scan the global 
structure of phase space in the spacetime (Al). This 
operation is realized by calculating thousands of orbits in 
practice. First we fix 9 = 0, meanwhile we let r run from 
r m in = 5.77 to r max = 22.34 with a span of Ar = 0.1. 
Then, for each given r we can solve for r and find two 
roots r_ and r+ (r_ < and r + > 0) from the Eq. 
(A4). Finally we take f from f_ to r + with a sampling 



interval Af = 0.01. Once r and f are given initially, 9 
should be derived by the relation (A4). As is shown in 
Fig. 5, we have plotted all the starting points on the r — f 
plane (9 = n/2). According to different values of FLI, we 
classify orbits by some dynamical features. This is called 
as the description of the global structure of phase-space 
for the system. 

By a comparison of Fig. 5 with Fig. 1, it is clear that 
the global dynamical features they depicted are nearly 
compatible. As an emphasis, an advantage to use FLI be- 
comes more apparent because it applies to systems with 
an arbitrary number of dimensions. Without doubt, it 
is handy to study the global dynamics of the spinning 
compact binaries [10] by virtue of FLI. It is also easier 
to probe the variation of the dynamical characteristics as 
certain parameters of the system vary. 



IV. SUMMARY 

We compared three methods for computing LEs with 
coordinate gauge invariance in general relativity. The 
three methods are the geodesic deviation vector tech- 
nique (Ml), the two-nearby-orbits method with projec- 
tion operations and with coordinate time as the inte- 
gration independent variable (M2), and the two-nearby- 
orbits method without projection operations and with 
proper time as the integration independent variable 
(M3). The contributions of this work are as follows. 

It is unnecessary to adopt the projection operators for 
Ml and M3 from the theoretical point of view. Method 
M3 is the simplest method for calculating LEs in a rel- 
ativistic gravitational system in most cases. As another 
contribution, we extended FLI to a coordinate invariant 
form in relativistic dynamics, and proposed its algorithm 
with the two-nearby-trajectories method. This indicator 
can rapidly, reliably and accurately distinguish between 
ordered and chaotic motions in relativistic astrophysics. 
Only when the initial separation and renormalization 
within a reasonable amount of time span are chosen ap- 
propriately, is this FLI nearly consistent with the FLI 
computed with the geodesic deviation vector technique. 
However, the former is rather simpler to use than the 
latter. As a characteristic, our FLI grows exponentially 
with time in the chaotic case, while it grows algebraically 
in the case of quasi-periodic trajectories. Evaluating the 
different behaviors of this FLI, we successfully explored 
the global dynamics of phase space where regions of chaos 
and order are clearly identified. 

It should be pointed out that a main advantage of the 
LE and FLI with two-nearby-trajectories is their easier 
application in treating complicated relativistic gravita- 
tional systems with high degrees of freedom, such as 
multi-body problems and spinning compact binaries. In 
the future work, we shall employ an appropriate numer- 
ical integrator and the FLI to find out what happens in 
the spinning compact binaries. The global dynamics of 
the binary systems will be explored, and the transition 
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from regular motion to chaos will be considered as some 
dynamical parameters of the system vary. 



Appendix A: Core-shell system and the equations of 
motion 

A Weyl spacetime including a non-rotating black hole 
surrounded by an axially symmetric shell in Schwarzchild 
coordinates (t, r, 9, <f>) is expressed as [6] 

ds 2 = -{l--)e p dt 2 + eQ- p [(l--)- 1 dr 2 
r r 



We now have to do the tedious derivation of the geodesic 
deviation equations. The result is too complex and will 
not be written here. If we take a' as the derivative of 
a with respect to coordinate time t, the above geodesic 
equations are readjusted as follows 



(A7) 
(A8) 



r 


r'di 


¥ 


~ Jdt' 


e 


9' di 


i 2 


i di' 



dt 
di 



+r 2 d9 2 } + e- F r 2 sin 2 9d<t> 



(Al) 



-Er(r-2)-^(fr' + ^e') 
-2E(r-2)- 2 e- p r'. 



where P and Q are the functions of r and 9 only. This 
system obviously has two integrals, the energy E and the 
angular momentum L in the forms 



i = Er(r-2)- 1 e- p , 
6 = Le p r- 2 sin" 2 9. 



(A2) 
(A3) 



In addition, the 4-velocity of a particle always satisfies 
the constraint 



U a U a = -(l-l) e p i 2 + eQ- p [{l-- r )-H 2 



+r 2 9 2 ]+e- p r 2 sin 2 94> 2 



-1. 



(A4) 



If a fourth integral holds, the system is integrable. Other- 
wise, it is possible to yield chaos. However, it is not easy 
to deduce whether the fourth integral exists or not. To 
study the dynamics of a geodesic particle in this system 
further, we need the following geodesic equations 



"(r 



1 1 dP dQ 

~~Y) + 2 { '~dr ~ ~dr~ 



)]f 2 



,dQ 8P ■ 



r ' i)Q nP )]9 2 +fl (r,9), (A5) 



dr 



2r{r-2) y d9 d9 ' 



2 dQ_dP ) . d 
r dr dr 



(A6) 



where 

Mr, 6) = 

Mr, 6) = 



1 



-)e 



2P- 



+ {!--)— ]t 2 

r or 

BP ■ 



+-(r - 2)e~ Q sir 
dP 

■sin9- ( — — - sin0 + 2cos0)0 2 . 

Ou 



Here r and 9 in Eqs. (A5) and (A6) should be changed 
into r — r'i and 9 = 9'i respectively, while i and <fi remain 
the original expressions of Eqs. (A2) and (A3) because 
they do not contain coordinate time t explicitly. 

Given two metric functions P and Q, the space- 
time (Al) is determined at once. Now we reexamine 
the full relativistic core-shell configuration involving a 
Schwarzschild black hole plus a pure octopole shell stud- 
ied by Vieira & Letelier [6] . In this model, we have 



P(u,v) 
Q(u,v) 



where u = r 
eter. 



-Ouv(5u 2 - 3)(5w 2 - 3), 
-lov[5(3u 2 - - v 2 ) - 4} 

+ ±- Q 2 [-2^\l-v 2 ) 

■{hv 2 + 2v - l){hv 2 - 2v - I) 
+15u 4 (l - w 2 )(65« 4 - AQv 2 + 3) 
-3?i 2 (l - v 2 )(25v 2 - 3)(5v 2 - 3) 
-v 2 {2Hv 4 -45v 2 + 27)], 

v = cos 9 and O is an octopolar param- 
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FIG. 1: Poincare surface of section at the plane = 5 (0 < 0) 
for the full relativistic core-shell configuration consisting of a 
Schwarzschild black hole and a purely octopolar shell with 
parameters E = 0.9679, L = 3.8 and O = 7.012 x 10" 7 . 




FIG. 2: LEs for a chaotic orbit with the initial conditions r — 9 and r = 0.025. The three methods, Ml, M2 and M3, give 
almost the same final value of LEs. 



11 




1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

iog 10 x iog 10 t 



FIG. 3: The relation of FLIs with proper time by use of the index (23) for the regular orbit M in the left (colored light gray). 
The black dots correspond to the result obtained from Eq. (22). The right is the same as the left but for the chaotic orbit N 
in Fig. 1. 




FIG. 4: Left: the evolution of the FLIs with the initial values of r ranging from the interval [5.77, 22.34] when the indicator 
(23) is used and r = is fixed at the initial time. Each orbit is integrated numerically till r reaches 10 5 . Two types of orbits 
consisting of regular and chaotic orbits are obtained according to distinct values of FLIs. Right: the same as the left, but r — 9 
is fixed first, and f takes values from the interval [—0.1428, 0.1428]. 
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FIG. 5: Regions of different values of the FLIs on the plane 
6 = ^. Initial conditions are colored white if their FLIs<6, 
and black if FLIs>6, which are corresponded to the regular 
region I, and chaotic region II, respectively. Here, all inter- 
section points of an orbit with the plane should be regarded 
to have the same final value of the FLI. 



